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Abstract Repaid global climate changes in 
temperature and rainfall influence the species 
distribution and diversity patterns. Chinse skink is a 
common species with large population and widely 
distribution in China. To access potential effect of 
climate changes on the unendangered species, we used 
the maximum-entropy modeling (MaxEnt) method to 
estimate the current and future potential distributions 
of Chinese Skink. Predictions were based on two 
periods (2050 and 2070), three general circulation 
models (GCMs: BCC-CSMI-1, HadGEM2-ES, MIROCS), 
four representative concentration pathways (RCP: 
2.6, 4.5, 6.0 and 8.0) and 28 environmental variables 
including topography, human impact, bio-climate and 
habitat. We found that the model were better fit with 
high values in AUC, KAPPA and TSS. The jackknife 
tests showed that variables of BIO9, BIO14, BIO15, HFI 
and GDP were relatively higher contributions to the 
model. Although the size of suitable areas for skink 
have less effect by future climate change under full and 
mull dispersal hypothesis, we should still focuse on the 
effect of human impact and climate changes on the 
protection and management for Chinese skink due to 
the variables uncertainty. 
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1. Introduction 


Many evidences indicated that rapid global climate change is 
already undergoing. Global average temperature has increased 
by 0.85°C from 1880 to 2012, and this trend is likely to continue 
for decades (Stocker et al., 2013). According to the most recent 
Intergovernmental Panel on Climate Change Fifth Assessment 
Report (IPCC AR5), the global average temperature may 
increase further by a minimum of 0.3-1.7°C (RCP 2.6) to a 
maximum of 2.6-4.8 °C (RCP 8.5) by the end of this century. 
Continuing global warming forced the species to response the 
climate change by expansion or contraction their native range, 
thus affected the global diversity by some biological events, such 
as habitat destruction, fragmentation and biological invasion 
(Pearson et al, 2007; Lyu and Sun, 2014; Penman et al., 2010). 

To investigate the changes of species distributions and 
diversity patterns under climate change, many species 
distribution models (SDMs) were developed and applied to 
detect the relationships between species distribution and 
environmental variables(Elith et al., 2006; Peterson, 2006). 
Although recent studies pointed out some drawbacks of 
SDMs, including neglecting competitive interactions, species 
plasticity, adaptation and time-lag (Davis et al., 1998, Hannah 
et al., 2002; Pearson and Dawson, 2003), SDMs could still be 
regarded as the useful models for predicting species potential 
distribution and extinction risk, assessing reserve designs, and 
evaluating conservation priorities by incorporating meta- 
population demography and landscape interactions, species 
life history traits, biotic interactions, species dispersal ability, 
species evolution and adaptation, and human activities 
and disturbances into the models (Keith et al., 2008; Preston 
et al., 2008; Li et al., 2014; Luo et al., 2014). 

Comparing to the mammals and birds, the poor dispersal 
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ability and autonomic thermoregulation lead to the reptiles 
are more sensitive to the climate change. Previous studies on 
reptiles were focused on the potential distribution prediction 
on some endangered and invasive species. However, the 
species with a stable population and wide distribution is less 
studied to our knowledge. Chinese skink (Eumeces|=Plestiodon] 
chinensis) is medium-large sized oviparous scincid lizard, widely 
distributing in East Asia including China, Korea and Vietnam 
(Lin et al., 2006). Based on the previous ecology studies, the main 
habitat for Chinese skink is lowland area including farmland, 
forest, and roadside (the altitude ranging from 10 to 1030 m), so 
the population size of this species is stable because of the strong 
environment adaptiveness (Pan et al., 2005). 

In this study, we selected the Chinese Skink as the model 
species using MAXENT model to simulate the potential 
distribution under different climate change scenarios for 
2050 and 2070. The aims of this study are: 1) to determine the 
current and future potential distribution of E. chinensis; 2) to 
find out the most important environmental factors in model 
predicting; 3) to compare the future distribution pattern of E. 
chinensis to the species with small population size and narrow 
distribution under future climate scenarios; 4) thus to provide 
recommendations for the protection and management of E. 
Chinensis in China. 


2. Methods and Materials 


2.1. Occurrence Data and Environmental Variables We 
collected species occurrence data from 3 sources: 1) network 
databases, including Global Biodiversity Information Facility 
(GBIF: http://www.gbif.org/), Specimen Resources Sharing 
Platform for Education (SRSPE: http://mnh.scu.edu.cn/main. 
aspx) and China Animal Scientific Database (CASD: http:// 
www.zoology.csdb.cn/page/index.vpage); 2) the records in the 
herpetological museum of Chengdu institute of biology (CIB), 
Chinese academy sciences (CAS); 3) the information on the 
published literatures (all the data accessed in December 2014). 
To ensure that only one record in each grid cell, we removed 
duplicate records in the same 1 km” grid cell. Totally 79 grid 
localities were kept for the following procedures. 

We collected 28 environmental variables and divided them 
into 4 categories: topography, human impact, bio-climate and 
habitat. The codes, source and resolution of each variable were 
shown in Table SI. The correlation between the candidate 
environmental variables was analyzed by the band collection 
statistics function in ArcGIS 10.2. We removed highly correlated 
variable (r 5 0.7) to minimize the impact of multicollinearity 
and over-fitting of the model (Burnham and Anderson, 2002). 
Finally, BIO3, BIO 5, BIO9, BIO14, BIO15, ELE, SLP, ASP, LAC, 
NDVI, POP,GDP, HFI were kept in the next model construction. 
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2.2. Global Change Scenarios For projecting future climate 
scenarios, we selected three internationally recognized general 
circulation models (GCMs): BCC-CSM1-1, HadGEM2-ES, 
MIROCS based on the Coupled Model Intercomparison Project 
Phase 5 (CMIP5: http://cmip-pcmdillnl.gov/cmip5). According 
to IPCC AR5 on greenhouse gas scenarios, we selected 
representative concentration pathways (RCP: 2.6, 4.5, 6.0 and 
8.0) in two time periods (2050 and 2070) (Moss et al., 2008). To 
reduce the uncertainty and illustrate future climate conditions, 
we averaged three GCMs for each RCP scenario to produce 
a total of 8 future models (four for the 2050s and four for the 
2070s). Because we could not projected future scenarios to the 
factors of habitat (AI, LAC and NDVI), human impact (POP, 
GDP and HFI) and topography (ELEV, SLP, and ASP) (Table 
S1), we assumed these factors were constant, and used them in 
these future predictions (Luo et al, 2014; Thuiller et al., 2006). 


2.3. SDM construction and evaluation We used MAXENT 
3.33 to prediction the Chinese skink distribution. The settings 
for the this software were as following: 1) regularization 
multiplier, 1; 2) maximum iterations, 500; 3) convergence 
threshold, 10°; 4) maximum number of background points, 
10 000; 5) random test percentage, 20%; 6) replicates, five; We 
selected the ‘auto features’ dependent on the number of presence 
records to reduce over-fitting (Phillips et al., 2006). We selected 
a logistic output format to keep the environmental suitability 
values ranging from 0 to 1, and carried out jackknife analyses 
of the regularized gain with training data to examine the 
importance of individual predictors. 

To evaluate the accuracy of each model, we firstly selected 
the AUC the area under the receiver operating characteristic 
(ROC) curve. AUC values ranged from 0 to 1, and models 
with values above 0.75 are considered potentially useful 
(Fielding and Bell, 1997). To test the reliability of the accuracy 
assessment, we then calculated Cohen’s kappa and true skills 
(TSS), which are both ranged from-1 to 1. Values 1 indicated a 
perfect performance, while values O indicated a performance 
no different to random (Burnham and Anderson, 2002). We 
used the presence.absence.accuracy function of the R package 
Presence Absence to calculate the kappa and TSS values (Freeman 
and Moisen, 2008). 


2.4. Impacts of Climate Change We first selected the 
thresholds by the method of minimum training presence and 
transferred raw outputs to presence absence maps (Pearson 
et al., 2007). Then two spread assumptions were used to 
calculate the range shift under climate change: 1) null spread 
(no spreadability of Skink); 2) full spread (unlimited ability to 
spread). Under the assumption of null spread, only the overlap 
habitat between current and future ranges was considered 


suitable for Skink. Under the full spread assumption, the Skink 
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populations could reach all new potential habitat ranges. We 
finally quantified the predicted current range (CR), potential 
range loss (RL, current suitable areas projected to be lost), and 
potential range gain (RG, current unsuitable areas projected to 
become suitable) by summing them by pixel, and calculated 
the RL and RG rates by dividing each by CR. Furthermore, 
we estimated the percentages of range change (RC) and range 
turnover (RT) using the following equations (Hu and Jiang, 
2011): 

RC = 100x (RG-RL)/CR 

RT = 100x (RG + RL)/(CR + RG) 


3. Results 


3.1. Model Performance and Variable Importance Values 
of kappa and TSS were both larger than 0.7, and the AUC 
values were larger than 0.9, which indicated that the SDMs 
were better fit for the Chinese Skink. And the five replicates 
for cross-validation were stable indicated that the model was 
robust (Table 1). The jackknife test showed that the training 
gain with all variables is 2.57 (n = 5), and them for each variable 
ranged from 0.0264 to 1.6125. Variables of BIO9, BIO14, BIO15, 
HFI and GDP achieved the highest gains when used in isolation 
indicated that they had relatively higher contributions to the 
model (training gain with only the variable > 1.10). However, 
other variables including ELEV, ASP, SLP, BIO3, BIO5, LAC 
and NDVI had limited contributions to the model (training gain 
with only the variable < 1.10) (Figure 1). Variable contributions 
analysis showed that BIO14 had the highest contribution (52.0%) 
to the model followed by HFI with a 27.7 % contribution, but 
POP and BIO15 had the lowest contribution (0.2%). 


3.2. Current and Future Potential Distributions and 
Range Shifts under Climate Change Based on the current 
suitability map, Chinese Skink is predicted with high habitat 
suitability in the southeast China with redder color, where is in 
the lower altitude region with more rainfall and higher mean 
temperature (Figure 2). According the threshold method of 10 
percentile training presence, the average values of the cutoff 
point is 0.156 (n = 5). With the threshold at 0.156, the potential 
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Figure 1 Results of jackknife test of relative importance of pre- 
dictor variables for Chinese Skink. 
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Figure 2 Mean predicted probability of occurrence (suitability) 
of Chinese Skink under current situation in China. The color rep- 
resents the suitability, from low (blue) to high (red). 


suitable area size for Chinese Skink is 1605446 km” at present, 
only accounts for 11.6% of the study area. 

To estimate the climate change impact, the potential 
distributions of Chinese Skink in 4 future climate scenarios 
were shown in Figure 3. Under the full dispersal hypothesis, the 
potential distribution in all four RCP scenarios is expected to 
increase by the 2050s (average range size = 1721108 km’, n = 4) 


Table 1 Accuracy measurements of predictive SDMs for Chinese Skink. AUC: area under relative operating characteristic curves; TSS: true 


skill statistic; Kappa: Cohen's kappa statistic, Rep.1-5 represent the five replicates for cross-validation. 


Rep.2 


Rep.1 Rep.3 Rep.4 Rep.5 
Accuracy Ensemble 
measurement Training Test Training Test Training Test Training Test Training Test 
kappa 0.774 0.800 0.783 0.794 0.724 0.798 0.631 0.803 0.756 0.803 0.847 
AUC 0.982 0.979 0.980 0.979 0.985 0.980 0.985 0.978 0.985 0.980 0.985 
0.766 0.783 0.799 0.703 0.741 0.523 0.803 0.713 0.803 0.847 


TSS 0.748 


-| Without variable 
With only variable ™ 
| With all variables = 
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Table 2 Areas of potential distributions and percentages of range loss (RL), range gain (RG), range change (RC) and range turnover (RT) of 


Chines Skink for Current, 2050, and 2070 under the climate change, based on the full dispersal and null dispersal hypothesis. 


Area (km’) RL (km) RG (km’) RC (%) RT (96) 

current 1605446 

Full dispersal hypothesis 

2050 rcp2.6 1692317 44541 131412 5.4196 10.13% 
rcp4.5 1767885 50514 212953 10.12% 14.49% 
rcp6.0 1682396 29164 106114 4.79% 7.90% 
rcp8.5 1741832 61520 197906 8.50% 14.39% 

2070 rcp2.6 1684729 46666 125949 4.94% 9.97% 
rcp4.5 1726442 50463 171459 7.54% 12.49% 
rcp6.0 1637441 74099 106094 1.99% 10.53% 
rcp8.5 1841842 51490 287886 14.72% 17.92% 

Null dispersal hypothesis 

2050 rep2.6 1560905 44541 0 -2.77% 2.77% 
rcp4.5 1554932 50514 0 -3.1596 3.1596 
rcp6.0 1576282 29164 0 -1.82% 1.82% 
rcp8.5 1543926 61520 0 -3.83% 3.83% 

2070 rcp2.6 1558780 46666 0 -2.9196 2.9196 
rcp4.5 1554983 50463 0 -3.1496 3.1496 
rcp6.0 1531347 74099 0 -4.62% 4.62% 
rcp8.5 1553956 51490 0 -3.2196 3.2196 


and further increase by the 2070s (average range size - 1722614 
km’, n = 4), respectively. However, the potential distribution 
under the null dispersal hypothesis showed an opposite 
tendency, all the potential distribution in 4 future RCP scenarios 
will decrease by the 2050s (average range size = 1559011 km’, n = 
4) and by the 2070s (average range size = 15497766 km’, n = 4) 
(Figure 3, Table 2). 

Across all scenarios in 2050s and 2070s, range loss (RL) 
under the full dispersal hypothesis were larger than that under 
the null hypothesis (Independent T test, t =7. 353, df = 14, P < 
0.01). Also from the values of range gain is 0 under the null 
dispersal hypothesis, because we assumed that the Skink could 
not expand their distribution by migration. Range gain (RG) is 
ranging from 106114 to 287886 km’, growing from 6.61%-17.93% 
to the current area, it indicated that the Skink will expand its 
suitable range under no limit on its dispersal. Range change is 
ranging from 1.99% to 14.72% under the full dispersal hypothesis, 
but ranging from -1.8296 to -4.62% under the null dispersal 
hypothesis, which indicated that the predicted potential suitable 
habitat for Skink present different results based on the species’ 
migrate ability. 


4. Discussion 


4.1. Sensitivity to Climate Change Many researches 


indicated that the global climate is undergoing rapid change, 
and is predicted to continue over the next century. The climate 
change on the earth is predicted to warm by 0.6 + 0.2 °C during 
the 20th century and by 14-58 °C in 21th century (Houghton 
et al, 2001). Warming of the earth may pose a threat to species 
by affecting population dynamics, distributions and the spatial 
structure of the suitable habitats (Bellard et al., 2012). Some 
studies predicted that changing climate will reduce the suitable 
habitat and narrow the distribution range, which will drive 
11% to 58% of vertebrate, invertebrate, and plant species to 
extinction by 2050 (Thomas et al., 2004). While some species are 
likely to benefit from the changes with extending ranges into 
currently unoccupied areas (Parmesan et al., 1999; Morueta- 
Holme et al., 2010). 

Most animal species may have at least limited dispersal 
ability and may not be able to colonize all of their climatically 
suitable area if other habitat requirements are not fulfilled 
(Levinsky et al., 2007). Full- and Null-constraint migration 
assumptions are the two extreme possibilities, and future range 
shifts will probably fall in between (Levinsky et al, 2007; Hu 
and Jiang, 2010; Luo et al., 2014). Chinese skink is a medium- 
large sized oviparous scincid lizard with ~134 mm snout- 
vent length (SVL), the population of this species is stable and 
wide distribution (Lin et al, 2006). Assume on the full spread 
hypothesis, our results predicted that the range loss (RL) of 
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Chinese skink is 2.89% (1.82% to 3.83%) by 2050s and 3.47% 
(2.91% to 4.62%) by 2070s, but range gain (RG) is 10.1% (6.60% 
to 13.26%) by 2050s and 10.77% (6.61% to 17.93%) by 2070s. 
We found that the range gain is much larger than the range 
loss for Chinese skink resulting from several climate change 
scenarios, which indicated that the current habitat of high 
occurrence probability is still be useful, and some current 
unavailable habitat could also be used in future. Even under 
the null dispersal hypothesis, Chinese skink were assumed that 
no migrate ability (RG = 0), the proportions of RL is relative 
small (2.89% by 2050s and 3.47% by 2070s). Most of the current 
suitable habitat is still with high occurrence probability by 
Chinese skink in future climate scenarios (RT < 15% in both 
hypotheses) (Table 2). Thus, depending on the population size, 
moment ability, suitable habitat, our modelling revealed that 
the suitable habitat of Chinese Skink was less affected by future 
climate change. 


4.2. The Importance of Environmental Variable 
Influencing the Future Distribution Previous researches in 
species distribution models to predict the potential distribution 
were focused on the bio-climate, topography and habitat 
factors (Pearson et al. 2007; Keith et al, 2008; Penman et al., 2010). 
Recently studies are paid more attention to the anthropogenic 
activity, because human activities represent the extensive, 
long-duration and persistent impacts on the environment 
that permanently alter ecosystem construction and ecology 
(Hu and Jiang, 2011, Li et al, 2014; Luo et al, 2014). As a result, 
we added the human impact into the MAXENT model. The 
Jackknife analysis in this study showed that the candidate 
factors in both bio-climate (BIO9, BIO14, BIO15) and human 
impact (HFI and GDP) groups influencing the distribution 
and patterns of range change (training gains > 10), however, 
factors in topography and habitat have limited effect on model 
prediction (training gain < 1.0) (Figure 1). Precipitation (BIO14 
and BIO15) and temperature (BIO9) represented the future 
climate changes in hydrothermal condition are benefit for 
Chinese Skink in prey and incubation (Ji and Zhang, 2000; Ji 
et al., 1995). 

Human impact (HFI and GDP) are other key factors 
contributed to the modelling prediction with training gains > 
LO (Figure 1). Previous ecology studies indicated that Chinese 
skink preferred the ecosystem changed by humans, such as 
farmland, grass and brushwood on the roadside (Lin et al., 
2006). And diets of Chinese Skink are mainly consisted to the 
invertebrates (annelid, molluscan and arthropod) covering 
more than 30 families (Lin and Ji, 1999). Human activity such as 
farming, harvest-picking and lumbering changed the habitat 
from coverage to open, which means that Chinese Skink 
could increase the detective probability on prey. Therefore, 
the increasing temperature, rainfall and human impact 
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will not only keep most of the existing suitable habitat for 
Chinese Skink, but also enlarge its range into some current 
unoccupancied region. Although this study does not necessarily 
provide accurate predictions of current and future distributions 
due to the small sample size and environmental data accuracy, 
it is a first model prediction on Chinese Skink under climate 
change and human activities in China and the results are still 
convinced that could be used to design the conservation plan 
for this species. 
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Appendix 


Table S1 Environmental Variables Used for Modeling of Chinese Skink (All the Variables were Compiled in WGS 1984). 


Variable groups Environmental Codes Data sources Resolution 
variables 
Annual Mean Temperature Biol Worldclim 1.4, http://www.worldclim.org/cmip5_30s ~1 km*(30-seconds) 
Mean Diurnal Pang (Meanofaionthly Bio2 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
(max temp - min temp)) 
Isothermality (BIO2/BIO7) (* 100) Bio3 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km?(30-seconds) 
Temperature Seasonality (standard Bio4 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
deviation *100) 
Max Temperature of Warmest Month Bio5 Worldclim 1.4, http://www.worldclim.org/cmip5_30s ~1 km*(30-seconds) 
Min Temperature of Coldest Month Bio6 Worldclim 1.4,http://www.worldclim.org/cmip5. 30s ~1 km”(30-seconds) 
Temperature Annual Range (BIO5-BIO6) Bio7 Worldclim 1.4, http://www.worldclim.org/cmip5_30s ~1 km*(30-seconds) 
Mean Temperature of Wettest Quarter Bio8 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
Mean Temperature of Driest Quarter Bio9 Worldclim 1.4,http://www.worldclim.org/cmip5_30s ~1 km’ (30-seconds) 
Bioclimate 
Mean Temperature of Warmest Quarter Biol0  Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km*(30-seconds) 
Mean Temperature of Coldest Quarter Bioll Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
Annual Precipitation Biol2 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km”(30-seconds) 
Precipitation of Wettest Month Biol3 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
Precipitation of Driest Month Biol4 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km”(30-seconds) 
Precipitation Ses rat (Craii Biol5. 8 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km(30-seconds) 
Variation) 
Precipitation of Wettest Quarter Biol6 9 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km”(30-seconds) 
Precipitation of Driest Quarter Biol7 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
Precipitation of Warmest Quarter Biol8 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km”(30-seconds) 
Precipitation of Coldest Quarter Biol9 9 Worldclim 1.4, http://www.worldclim.org/cmip5. 30s ~1 km’ (30-seconds) 
Elevation (m) ELEV Global digital elevation model, http://srtm.csi.cgiar.org/~1 km” (30-seconds) 
7 Calculated based on ELEV by the slope 3 2 
Topography Slope (°) SEP! function im AreCIO 102 1 km (30-seconds) 
` Calculated based on ELEV by the aspect 2 2 
Aspect () HOP function in ArcGIS 10.2 1km. (30-seconds) 
CGIAR-CSI Global-Aridity and Global-PET Database, 
Aridity index Al http://www.cgiar-csi.org/ ~1 km”(30-seconds) 
data/global-aridity-and-pet-database 
; ; Global Landcover 2000, http://ies.jrc.ec.europa.eu/ 2 
Habitat Habitat Landcover type LAC Bababa ena 2000 1 km (30-seconds) 
Normalized difference vegetation index NDVI Thematic Database of Human-Earth System, http:// ~1 km” (30-seconds) 
www.data.ac.cn/ 
Human population (individuals/km’) POP ‘Thematic Database of Human-Barth System, http:// ~1 km” (30-seconds) 


Human impact 


Gross domestic product (104 Chinese Yuan/ GD 


km’) 


Human footprint index 


P 


HFI 


www.data.ac.cn/ 


Thematic Database of Human-Earth System, http:// 
www.data.ac.cn/ 


Last of the Wild Data, http://sedac.ciesin.columbia.edu/ 


wildareas/ 


~1 km’ (30-seconds) 


~1km/’(30-seconds) 


